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THE  STABILITY  OF  COMPLICATED  ENERGY  SYSTEMS 

by 

V.  F.  Kurov 


The  original  equations  for  the  disturbed  motion  of  complicated  automatically  con¬ 
trolled  systems  can  be  represented  in  the  following  matrix  form: 


pX  =  AX\  X 


(1) 


where  the  differential  equation  (1)  is  the  first  approximation  of  the  mathematical  model 
of  the  system 


n 

P*t  =  /l  ^  Ol/-*/  +  9l  (•*! . •*,) 


/’•*»=/.  {•*! . ■«»)  = 


=  ^  «./■*/ +  9. 


(X|,  . .  .  ,  Xf^] . 


(2) 


Thus,  when  one  takes  into  account  the  basic  determining  factors,  the  solution  of  the 
question  regarding  stability  "in  the  small"  (static  stability)  reduces  to  analysis  of  the 


FTD-HT-23-888-67 


1 


system  (1)  of  linear  differential  equations  describing  the  electromechanical  and  electro¬ 
magnetic  transfer  processes.  As  usual,  the  characteristic  equation  of  the  system  in 
question  has  a  rather  high  order  and  can  be  represented  in  the  form 


S(/>)  =  Vfl,p""  =  0,  (3) 

o 

where  .... 

The  difficulty  in  analyzing  such  a  system  is  due  both  to  the  calculation  of  the  co¬ 
efficients  of  the  characteristic  equation  and  to  the  application  of  the  stability  criteria 
themselves,  which  are  quite  laborious  for  equations  of  high  degree  and  necessitate  spend¬ 
ing  a  considerable  amount  of  time  for  their  verification.  The  application  of  electronic 
digital  conq)uting  machines  requires  the  construction  of  optimal  algorithms  that  will  en¬ 
able  one  to  answer  the  above-mentioned  questions  in  the  least  amount  of  time  and,  in 
some  cases,  with  the  least  demands  on  the  memory  of  the  machine.  Here,  it  is  desir¬ 
able  to  use  ready-made  standard  programs  of  linear  algebra. 

In  the  present  article,  we  shall  consider  the  possible  methods  of  constructing  re¬ 
gions  of  stability  of  a  linear  automatic-control  system  in  the  space  of  the  control  param¬ 
eters  fay  using  electronic  digital  computing  macl^es  and  also  methods  for  obtaining  the 
cbe^cients  ctf  the  characteristic  equation. 

The  coefficients  of  the  characteristic  equation  (3)  are  independent  of  the  time  but 
they  depend  on  the  set  of  control  parameters  k^,  . . . ,  k^,  which  vary  in  some  closed 

bounded  region  K,  The  investigation  problem  consists  in  finding  a  closed  region  GQ  1( 
in  which  the  roots  of  the  characteristic  equation  will  have  only  negative  real  parts. 

In  searching  for  such  a  re^^on  (the  region  of  stability),  we  usually  use  the  methods 
of  a  D-partition  Il,2].  If  the  number  of  control  parameters  does  not  exceed  2,  then  the 
boundary  of  the  region  of  stability  can  be  obtain^  in  parametric  form  in  a  comparatively 
simple  manner  [3]. 

Althouc^  it  is  possible  in  theory  to  apply  the  method  of  a  D-p«irtition  for  a  larger 
number  of  parameters  [4],  the  calculations  in  such  a  case  involve  considerable  difficul¬ 
ties,  especially  if  the  degree  of  equation  (3)  is  high.  The  analysis  of  the  results  obtained 
is  considerably  complicated.  Althou^  use  of  an  electronic  digital  computing  machine 
alleviates  the  computational  difficulties,  the  necessity  of  changing  the  frequency  over  a 
wide  range,  the  necessity  of  choosing  the  optimal  step,  etc.  render  uneconomical  the  use 
of  such  a  machine. 

Therefore,  the  tendency  to  use  directly  the  criteria  of  Routh,  Hurwitz  and  Yu.  I. 
Naymark  [1,2]  directly  to  analyze  the  stability  and  to  find  boundaries  of  stability  in  the 
region  of  the  control  parameters  [5,6]  is  completely  justified. 

The  methods  of  constructing  regions  of  stability  that  are  considered  in  the  article 
are  based  on  application  of  algebraic  criteria  for  stability.  Therefore,  as  a  preliminary, 
let  us  look  at  ^e  possible  methods  of  calculating  the  coefficients  of  the  characteristic 
equation. 


k^),  where,  in  turn,  k^,  . . . ,  k^  are  the  control  parameters. 
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Methods  of  ObtaininR  the  Coefficients  of  the  Characteristic  Ec 


The  characteristic  equation  <  disturbed  motion  of  a  complicated  energy  system  can 
usually  be  represented  in  the  form  of  the  determinant  of  an  ^  th-degree  matrix  nolynomial 
(7.81. 

5i(P)  =  M./>'  +  ^p'"’  + +  =  (4) 

where  A^.  ,  A^  are  square  matrices  of  order  m. 

To  set  up  a  program  for  calculating  the  stability  and  the  electro-mechanical  trans¬ 
fer  processes  of  a  complicated  energy  system  on  a  electronic  digital  computing  machine, 
it  is  necessary  to  represent  this  determinant  in  the  form  of  a  scalar  polynomial.  To  this 
end,  we  can  use  algorithm  3  of  two  basic  forms.  The  first  of  these  includes  algorithms 
based  on  linearized  equations  of  a  disturbed  motion,  which  are,  as  a  preliminary,  re¬ 
duced  to  a  form  solved  for  the  derivatives.  Then,  the  coefficients  of  the  characteristic 
equation  can  be  found  by  e^anding  the  determinant  of  the  "secular"  equation. 

The  second  form  of  algorithms,  which  do  not  involve  the  necessity  of  reducing  the 
original  equations  to  normal  form,  was  developed  by  L.  V.  Tsukemik  [9-11]. 

The  methods  of  representing  the  determinant  (4)  in  the  form  of  a  polynomial  include 
the  direct  e3q>ansion  of  it  by  means  of  an  expansion  in  terms  of  the  elements  of  any  row 
or  column  and  also  the  reduction  of  the  determinant  (4)  to  a  triangular  form.  However, 
these  methods  require  either  a  large  number  of  operations  or  a  complicated  program.! 

Let  us  look  at  some  of  the  possible  cases  of  the  obtaining  of  the  coefficients  of  the 
characteristic  equation. 

1.  Transformation  of  Equation  (4)  to  the  Form  |B  -  XE I  =  0. 

To  obtain  the  characteristic  polynomial,  we  can  use  an  artificial  device  based  on 
the  theorems  of  matrix  algebra. 

An  mth-order  determinant  whose  elements  are  rth-degree  polyiiomials  can  always 
be  transformed  into  a  determinant  of  order  mr  whose  elements  are  linear  in  p  provided 

the  matrix  of  the  coefficients  of  p^  is  nonsingular. 

Siqppose  that  the  given  equation  is  of  the  form  ^4).  Then,  its  roots,  as  one  can 
verify  directly,  are  equal  to  the  characteristic  number  of  the  matrix 


• 

To 

1 

-Ao'Ar-t 

K, 

1  o 

1 

E 

0 

0 

0 

0 

0 

E. 

0 

(5) 

where  A  is  the  inverse  of  the  matrix  A^  and  E  is  the  unit  matrix  of  mth  order. 
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The  characteristic  equation  of  the  matrix  (5)  can  be  found  by  means  of  Danilevskiy^s 
method.  Because  of  the  special  form  of  the  matrix  (5),  the  calculating  process  is  car¬ 
ried  out  considerably  more  rapidly  than  with  other  methods  [12]. 

The  chariEicteristic  equation  of  the  matrix  (5)  can  be  represented  in  the  form 


n 


If  we  rqilaoe  X  with  p  in  equation  (5a),  we  obtain  the  characteristic  equation  of  the  sys¬ 
tem  in  question. 

However,  in  a  number  of  cases,  the  matrix  in  equation  (4)  may  not  have  an  in¬ 
verse.  If  we  replace  p  with  in  equation  (4),  we  obtain  the  equation 

I  -f-  Af  - 1  Pi  ^  -f-  . . .  A©  1  =  0,  (6) 

the  roots  of  which  are  obviously  equal  to  the  reciprocals  of  the  characteristic  numbers 
of  the  matrix 


—  Ar  Af  _  1  —  Ar  Af  -  2  ...  —  Ar  *  A|  —  Ar  *  A^ 

E  0  ...  0  0 


where  A’|;  is  the  inverse  of  the  matrix  A^,,  which  we  can  always  make  nonsingular  by 
choosing  suitable  values  for  the  control  parameters. 

The  characteristic  equation  of  the  matrix  (7)  can  also  be  represented  in  the  form 
of  equation  (5a).  Then,  the  characteristic  equation  of  the  automatic-control  system  is 
obtained  from  equation  (5a)  by  replacing  X  with  l/p. 

2.  Application  of  the  Method  of  Interpolation  to  Equation  (4). 


Since  the  function  S^(p)  in  equation  (4)  is  a  polynomial  of  degree  not  exceeding  n, 
the  interpolational  polynomial  S(p)  must  be  identically  equal  to  S^(p).  This  follows  from 
the  theorem  on  the  uniqueness  of  the  interpolational  polynomial. 

Expansion  of  the  mth-order  determinant  (4)  yields  a  polynomial  of  degree  n  ^  mr 
in  p.  This  polynomial  contains  n  +  1  coefficients  of  the  different  powers.  These  coef¬ 
ficients  can  be  determined  if  we  know  the  values  of  S^(p)  for  n  1  values  of  p.  To  do 

this,  we  can  use  the  interpolation  formula  or  we  can  solve  a  system  of  linear  equations 
for  the  coefficients.  If  the  difference  between  two  successive  values  of  p  is  constant. 


4 


then,  to  obtain  a  polynomial  expression,  we  can  use  the  first  interpolation  formula  of 
Newton  [12]. 


S  {p)  =  S  (»)  + -^  A  S  (i)  +  4*  S  (»)  +  . . . + 

+  *(*-»)  •••(*->»  + I)  5 

ii! 


(8) 


where  x  =  (p  -  i  )/h  and  the  differences  A*  5  (3)  can  be  taken  from  the  following  Table: 


S(») 

S{i  +  h) 

A5(«) 

A»S(«) 

5  (2  +  2A) 

A5(«4-A) 

A‘S(S  +  A) 

A5(«4-2A) 

S  (8  +  3A) 

• 

5  (8  4“  ^^) 

• 

where  AS(6)-5(6+A) -5(A).  A*S(A)-A5(6+^) -iiS(6),  etc. 

The  kth-order  difference  of  a  kth-degree  polynomial  is  equal  to  zero  [sic],  which 
enables  us  to  determine  immediately  the  order  of  the  characteristic  equation. 

Formula  (8)  does  not  immediately  yield  a  polynomial  form  since  each  term  in  it  is 
a  polynomial  in  p.  Let  us  consider  the  interpolational  polynomial  (8)  arranged  in  increas¬ 
ing  powers  of  p  -  A,  so  that  the  coefficients  can  be  calculated  immediately  with  the  aid  of 
Table  1.  ^  To  do  this,  we  transform  it  to  the  form 


H 


=  =  +  (p-«)+ Mp -»)•+..;  + 


+  (p  -  »)", 


(9) 


where 


*.  =  -jr  5  (*>•  (*  =  1. 2 . «).  p,  =  s (»).  (10) 


For  convenience  in  making  the  calculations,  we  need  to  set  A  =0  and  h  =  1,  which  leads 
to  a  certain  simplification  of  the  expressions  (9)  and  (10). 
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If  6  is  chosen  different  from  zero,  the  shift  to  the  degree  p  instead  of  (p  -  ft)  is 
easily  carried  out  with  the  aid  of  Homer’s  scheme  or  in  accordance  with  the  scheme 
shown  below: 

=  ««  -  1  fl  -  1  ft*  -  1  4-  ««  -  1  « 


Ufl  -  1  =  *11  ftl  +  «12  ftj  ®  +  •  •  •  +  «1  n  * 

=  «00^f +  *f  1  +  •  •  •  +  *•«  ^n®" 


Or,  in  the  general  case, 


**(*  +  »)  ft*  + 1  ft  =»  0, 1, 2, ...» a. 


Here, 


«•/=(- IK  «,/=0.  i>j\ 


k^Tti  sAii. 

i<A  (U«0.1 . n). 

The  values  of  the  coefficients  for  equations  of  the  first  sixteen  orders  are  shown  in 
Table  2. 

In  the  case  of  constant  differences  between  successive  values  of  p,  Newton's  first 
interpolation  formula  becomes  unsuitable,  in  this  case,  we  can  use  Lagrange's  formula 
or  the  method  of  undetermined  coefficients. 

3.  The  Application  of  the  Method  of  Undetermined  Coefficients. 

Let  us  look  at  an  expression  of  the  form 

S(p)  =  a„  +  On-iP -f- •  •  •  +  <*iP"  *+^oP"- 

If  we  choose  n  +  1  values  of  p,  then  to  determine  the  coefficients  we  need  to  solve 

r  +  1,  linear  equations  (or  n  linear  equations  in  the  case  in  which  one  of  the  values  of 
p  is  zero,  in  which  case  S(0)  =  a^) 


a,  +  +  a._»pi  +  ...  +  a,pf“'  +a,/>"  =  S  (/»,), 

(/=:t.2,....«  +  l), 
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or,  in  matrix  form, 


where 


P  •  a-A 


(13a) 


1  Pt  Pt  ...  p’ 

•  Pi.  pi  •  Pi 

»  Pt  Pi  •••  pi 

a  = 

flu -2 

S  = 

S{p,) 

S{p,] 

1  Pm  pl  •••  pi 

flo 

sIpm) 

(14) 


The  system  of  equations  (13)  with  p.  =  S  +  ih  can  be  reduced  by  means  of  elementary 
transformations  [12]  to  the  equivalent  step  form 


P'-‘«  ‘6  =  4S(«).  *0  = 


» 


where  the  coefficients  b  determine  the  characteristic  equation  S  (p  —  i)*  =  0.  The 

i  s  0 

imaginary  axis  in  the  complex  plane  of  the  roots  of  this  equation  is  displaced  to  the 
right  (resp.  to  the  left)  for  positive  (resp.  negative)  8  with  respect  to  the  imaginary 
axis  of  the  complex  plane  of  the  roots  of  ^e  equation 


2a,/-'  =  0. 

i  s  0 


Here,  is  the  matrix  inverse  to  P*,  H  =  [h  h  ...  h*']  is  a  diagonal  matrix,  AS  (8)  = 

=  (A S  (8)*A*S  (S) . . .  A''S(B)}  and  A  =  (Aj  82 . . .  are  the  column  matrices  of  the  correspond¬ 
ing  quantities. 

One  cm  easily  see  that  finding  the  coefficients  of  the  characteristic  equation  by 
using  the  inteipolation  method  (10)  is,  by  its  nature,  a  more  rational  procedure  for 
solving  the  problem  than  the  method  of  undetermined  coefficients  in  the  case  of  equally 
space  values  of  p. 


To  solve  the  system  of  equations  (13a),  we  can  use  the  familiar  methods 

a«P  'S.  (13b) 

In  the  case  in  which  the  interval  between  adjacent  values  of  p  is  equal  to  1  and  p^  ~  0, 
the  system  of  equations  can  be  written  in  the  form 


(13c) 
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The  matrix  P  in  the  ejqpreeaion  (13c)  is  ind^endent  of  the  numerical  values  of  S(p  ).  We 

-1  ^ 

can  calculate  in  advance  the  inverse  matrix  P  ,  and  then  the  solution  of  ^e  system  (13c) 
reduces  to  multiplying  matrices  in  accordance  v^th  the  e^qpression  (13b). 

Figure  1  shows  the  number  of  operations  necessary  to  calculate  the  characteristic 
polynomial  different  methods  fcr  r  =  2,  depending  on  the  order  of  the  matrix  m. 


••Ip' 


Fig.  1.  The  dependence  of  the  number  of  com¬ 
putational  operations  on  the  order  of  the 
determinant. 

1)  Direct  ejqsansion;  2)  The  method  of  undeter¬ 
mined  coefficients;  3)  The  method  of  the  in¬ 
terpolation  formula;  4)  Tr;insformation  to  diag¬ 
onal  form;  5)  The  method  of  undetermined  co¬ 
efficients  with  a  constant  stq>  when  there  is  an 
inverse  matrix. 


As  one  can  see  from  Fig.  1,  direct  e^ansion  yields  the  smallest  number  of  opera¬ 
tions  if  the  determinant  is  of  fourth  order;  the  interpolational  method  and  transformation 
to  diagonal  form  are  preferable  for  hi^er  orders.  The  interpolational  method,  however, 
is  the  most  convenient  one  even  if  it  requires  somewhat  more  operations.  This  method 
enables  us,  by  the  nature  of  the  change  in  the  function  of  the  differences,  to  predict  the 
order  of  the  characteristic  equation.  The  values  of  the  function  can  be  calculated  fairly 
easily  with  the  aid  of  an  electronic  digital  computing  machine  with  the  minimal  number  of 
iiq)ut  data.  Althou^  the  method  of  undetermined  coefficients  requires  fewer  operations 
in  the  case  of  equally  spaced  values  of  p  than  does  the  method  of  transformation  to  diag¬ 
onal  form,  the  necessity  of  finding  the  inverse  of  a  matrix  of  high  order  makes  this  meth¬ 
od  less  desirable  when  one  is  using  an  electronic  digital  computing  machine. 
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It  is  convenient  to  employ  the  method  of  transformation  to  diagonal  form  when  one 
has  foimd  even  ’’inexact”  thou^  easily  visualized  regions  of  stability  I13J.  (By  ’’inexact” 
regions,  we  mean  regions  somewhat  smaller  than  the  usual  regions  of  stability  defined 
by  the  Routh-Hurtz  inequalities. ) 

Thus,  the  use  of  Newton’s  interpolation  method  is  extremely  convenient  when  one 
is  programming  an  electronic  digital  computing  machine  since  it  makes  it  possible  to 
use  the  standard  programs  of  linear  algebra.  In  addition,  this  method  compares  favor¬ 
ably  with  the  others  in  that  it  makes  fewer  demands  on  the  memory  of  the  machine  at 
the  cost  of  a  row-by-row  formation  of  the  matrix  of  the  transformation. 

Let  us  indicate  the  order  in  which  we  obtain  the  characteristic  polynomial  if  the 
control  parameters  k^ . k^  appear  linearly  in  equation  (3),  that  is,  if  the  control 

is  carried  out  with  respect  to  the  individual  parameters. 

If  we  set  k^,  . . . ,  k^  equal  to  zero  and  reduce  the  determinant  (4)  to  polynomial 

form,  we  obtain  the  coefficients  of  the  characteristic  equation,  which  are  independ¬ 
ent  of  K,  If,  in  the  expression  (4),  we  set  all  K  equal  to  zero  with  the  exception  of 
one,  we  can  obtain  the  coefficients  of  the  characteristic  equation,  which  depend  only 
on  ihe  given  parameter,  etc. 3 


FOOTNOTES 

(p,  3)  ^Furthermore,  to  reduce  the  determinant  to  triangular  form,  one  can  change  the 
order  of  the  characteristic  equation  at  the  cost  of  computational  errors. 

2 

(p,  5)  The  coefficients  P*  (d  can  be  calculated  from  the  following  recursion  formula: 


P*(i  +  i) 


P»-l(i)  — *PlKf) 

«  +  1  ’  (1) 


A. 


When  we  use  an  electronic  digital  computing  machine  to  calculate  the  coeffi¬ 
cients,  a  more  convenient  formula  is 


2  |P*-iwl 

Pm* +  !)-<-  ■■  - ;a-0,  i,2...  .  n;  *-1,2,. ...a. 


Here, 


Po  (0)  "  *  •  PO  (I)  ^  =  1 1  2, ... ,  A. 

3 

(p.  11)  If  the  coefficients  of  the  characteristic  equation  depend  on  constants  of  the  time 
of  the  regulators,  on  the  system,  or  on  other  parameters,  the  determination  of 
tlie  coefficients  of  the  characteristic  equation  as  a  function  of  these  parameters 
must  necessarily  be  carried  out  with  consideration  of  the  fact  that  the  degree  of 
the  characteristic  polynomial  will  change. 
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